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Abstract 

We study thermal Casimir and quantum non-retarded Lifshitz interactions between dielectrics in 
general geometries. We map the calculation of the classical partition function onto a determinant 
which we discretize and evaluate with the help of Cholesky factorization. The quantum partition 
function is treated by path integral quantization of a set of interacting dipoles and reduces to a 
product of determinants. We compare the approximations of pairwise additivity and proximity 
force with our numerical methods. We propose a "factorization approximation" which gives rather 
good numerical results in the geometries that we study. 

PACS numbers: 05.10.-a 41.20.Cv 
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I. INTRODUCTION 



Much is known about fluctuation-induced interactions between bodies. One distineruishes 
.wo pans to the ,„te..actio„ g The fl.t . ,ua„tun, in nature, and co^esponds 

to a multi-center generalization of London dispersion forces. The second, the so called 
thermal Casimir effect, has as its origin the thermal excitation of polarization modes in 
dielectrics, giving rise to temperature dependent forces. Together they form the multi-body 
generalization of the long-ranged part of the van der Waals interaction. An alternative, but 
equivalent vision comes from associating these forces to the energy and entropy of fluctuating 
electro dynamic flelds. A number of sophisticated theoretical techniques have been applied 
to the calculation of this interaction, but only the simplest of geometries are analytically 
tractable to exact solution, for instance planar surfaces 4|], spheres or cylinders. 

More complicated physical situations are generally studied with perturbation theory, or 
in exceptional cases with non-perturbative methods adapted to individual geometries j^. 
Such methods, that keep into account the totality of the interaction, are usually developed 
on the basis of Lifshitz theory for dielectrics. Analytic methods are typically based on 
an approximation around one of the few solvable problems. These range from the simple 
proximity force approximation, to a semi-classical interactions approximation [6], to the 
multipole expansion [7], to the more recent and sophisticated ray optics approximation ^. 
The main problem with these methods is they are limited when the system under study 
departs too far from the unperturbed system. Given the limitations of all analytic methods, 
which are in practice only capable of studying very regular geometries, numerical approaches 
have been developed. These methods have been limited both by storage space required and 
by long computational times, making it possible to study only two-dimensional systems, 
or three-dimensional systems translationally invariant in one direction [o], [lo|. It is only 
very recently that more powerful methods have been developed based on the numerical 
calculation of the Maxwell stress tensor Uj. The numerical approach we present in this 
paper is complementary to this last method, but it is based on somewhat different theoretical 
grounds. We will explain in the discussion how the two methods compare in detail. 

We note, as an aside, that so called "atomistic modeling" of materials, as performed with 
most molecular dynamics codes, with the aim of understanding micro-mechanical response, 
conformations of macromolecules or interfaces, often uses assumptions such as pairwise ad- 
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ditivity of interactions 12(] parameterized with phenomenological Lennard- Jones potentials. 
These methods neglect retardation, screening of the classical interaction by ions and multi- 
body effects. It is of clear interest to develop methods which will enable one to have a 
better quantitative understanding of such effects in both soft and hard condensed matter 
physics. More sophisticated (and even more costly) quantum simulation methods based on 
local densit y fu nctionals exsist, but are known to miss long-ranged dispersion interactions 



completely 



13|. 



In this paper we want to focus on important geometries which are clearly beyond per- 
turbative study. One simple example is the case of a tip near a structured surface, where 
the curvature of the surface makes most approximation methods ineffective, unless the di- 
electric contrast is very weak pjj]. This paper has as its principle aim the generalization 



of a recent paper 



15| which used direct diagonalization of a large matrix in order to calcu- 



late the free energy of fluctuating dielectrics in the classical (thermal Casimir) limit. The 
methods in [l^ were rather limited, only very small three dimensional systems could be 
studied. Quantum effects were neglected entirely. In this paper we use more sophisticated 
factorization techniques which allow us to treat finer discretizations of physical systems. We 
will, in addition, show how to introduce quantum mechanics into our formalism and use it 
to study the full, non-retarded interaction between two dielectric bodies. In order to com- 
pare with our numerical method, we analyze the performance of two simple approximation 
methods: pairwise additivity, and the proximity force approximation. We also introduce 
an approximation based on the factorization of the geometric and material properties that 
works surprisingly well for the geometries that we study. 

We begin by formulating the theory of classical fluctuating dielectrics. In this we prefer 
a formulation of the classical interaction in terms of the true microscopic fluctuating field, 
the polarization, rather than the electric field and electric displacement. We then show how 
to efficiently factorize the resulting quadratic forms and apply our formalism to calculate 
the free energy of interaction between a tip and an indentation. We then generalize our 
microscopic energy functional in order to consider the quantum fluctuations of a dielectric 
and evaluate the non-retarded interaction between a tip and a structured surface. 
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II. CLASSICAL FLUCTUATING DIELECTRICS 



We begin by evaluating the classical thermal interaction, since most of the technical 
difficulties of discretization and factorization are already present, before showing how to 
treat the quantum case. We start with the energy of a heterogeneous dielectric system 
written in terms of the polarization density p. We note that formulation in terms of the 
polarization field is also much more convenient for the generalization to scale dependent 
dielectric effects 17|, which are particularly important in water based systems. 

The energy g-ranged Coulomb 

energy for the induced charged density pi = —V- p and secondly a local contribution which 

depends on the local electric susceptibility Xo(r). 

/• V-p(r)V-p(r')^3 , f p'(r) ^3 

U„ = / i-V ^-^(Tr dV + / ^ \ \ . (1) 

J 87r|r-r'| J 2xo(r) ^ ^ 

We note that the dielectric constant e(r) = 1 + Xo(r); we use units where eo = 1. The 

partition function of the fluctuating dipoles 

Z = Jvp expi-pU,), (2) 

can be simplified by re- writing the Coulomb potential, 1 / Airr as an integral over a potential 
0. We then find the effective energy 



^P4 



One now performs the Gaussian integral over p, to find the partition function expressed as 
an integral over 0. To do this it is useful to integrate by parts replacing V- p with — p- V 0; 
boundary terms vanish in periodic systems or those in which fields decay at infinity. 

In order to perform numerical calculations, we derive a discretized theory. We discretize 
by placing scalar quantities such as (p on V = nodes of a cubic lattice. We choose a 
length scale such that the lattice spacing is unity. The lattice spacing should be much larger 
than the atomic scale, so that a formulation in terms of continuum dielectric properties is 
possible; the lattice should however be sufficiently fine to resolve features of physical interest, 
such as points or rough surfaces. Vector fields, such as V 0, are associated with the 3V links, 
e and xo are also associated with the links. We find 

Z{e) = n Xof I 'D<P exp (^(3 ^ |(9,0)^ j (4) 



where / is a link; the discretization of the derivative di evaluates the difference between 
variables on the corresponding nodes. The integral in eg. (jl]) is over all modes excluding 



the uniform mode, = const which has eigenvalue zero 15|. We write the quadratic form 
appearing in eq. (jlj) as (j)M(f)/2 with a symmetric matrix Af . The interesting, long ranged, 
part of the free energy of interaction comes from 

^=^log(det'(M)) (5) 

where the determinant again excludes the zero eigenvalue. We will see that the advantage 
of working with a matrix of the form eq. (jl]) rather than eq. ([1]) is that the quadratic form 
eq. (jl]) is sparse; a very large majority of its elements are zero, enabling the use of efficient 
evaluation strategies. While eq. (jl]) is also Gaussian, its quadratic form is dense. 

Some care has to be exercised to separate various contributions to the free energy eq. (jSj). 
In most situations of interest we want to find interfacial energies, which are subdominant 
compared with the bulk free energy of the full three dimensional system which we treat. 
The problem of separat ing bulk and surface contributions in general geometries is treated in 
a separate publication 15[. In the first physical example presented in this paper the volume 
of each phase changes as the tip is moved farther from the indentation, and the surface free 
energy is obtained after subtracting the bulk free energy. In the second system we analyze, 
the volume of each phase is constant so that variations in the interfacial free energies can 
be found by simple subtraction, as the bulk contribution cancels out completely. 



III. MATRIX FACTORIZATION 



In the simplest geometries, such as multiple parallel plates, the determinant is evaluated 
analytically from the eigenvalue equations for sets of plane waves j^. Since we wish to work 
in arbitrary geometries we evaluate the determinant eq. (j5]) using general matrix methods, 
with periodic boundary conditions to minimize edge effects. Despite the fact that M is 
sparse, and can be stored in memory which is proportional to the volume V = L^, standard 
(dense) matrix methods find the eigenvalues in a time which varies as Tg ~ V^'^, using 
a storage 5 ~ V"^. Both these scalings turn out to be prohibitive, and limit us to the 
study of systems no larger than linear dimension L = 25. We thus turn to methods which 
are adapted to large sparse matrices. For such matrices it is preferable to evaluate the 
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determinant without forming intermediate dense matrices, saving computer memory and 
accelerating calculations. We evaluate the determinant by firstly modifying the matrix in 
order to render it positive definite, so that we are not troubled by the zero eigenvalue. We 
do this by adding V to a single arbitrary diagonal element of the matrix; this can be shown 
to be equivalent to neglecting the zero mode of the determinant p^]. With the matrix, M 
now positive definite we write it as a product of Cholesky factors M = R^R where R is an 
upper triangular matrix. Because R is triangular the determinant is given by the product 
of the diagonal elements. From the factor R we find det M = (det i?)^. 

This factorization has a number of remarkable properties that make it far more powerful 
than diagonalization. A good choice for the ordering that is used for evaluating the Cholesky 
factors renders the method particularly interesting: Nested dissection numbers the nodes 
in a non-consecutive manner by recursively cutting a graph into equal pieces [isl. For this 
ordering of our matrix the Cholesky factor (in three dimensions) can be generated with a 
storage of 5* ~ V^^^, and with arithmetic effort Xg ~ V"^ [3]. We used the software package 
Taucs 20|] to perform the Cholesky factoring. We find that a system of ^ = 64^ can be 
factored on a 3.2GHz 64-bit Xeon workstation in approximately 300 seconds, using 3GB of 
memory. 



IV. INTERACTIONS BETWEEN A TIP AND INDENTATION 

We now consider the interactions between a tip and a surface indentation in a system 
in which the thermal contribution is dominant. Experimentally this can occur either with 
tips in water which are made with materials which have similar optical properties (so that 
the high frequency contributions to the Lifshitz energy cancel out), or in systems that are 
sufficiently separated such that the quantum interactions have died out requiring separations 
larger than hc/2kBT ~ 5/im at room temperature 

We evaluated the free energy of interaction for a system composed of a three dimensional 
rounded tip close to an indentation in a surface discretized to a lattice of dimension L = 80^, 
Figure [H We evaluated the free energies with two vertical displacements of / = 6 and 
/ = 20. The physically interesting case for optically matched systems has a contrast ratio 
of approximately = ei/e2 = 50 between the two media. However in order to study the 
evolution of the interactions with contrast we worked with a minimum ratio of = 1.05 to 
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a maximum of = 100. 




FIG. 1: Section of the tip-indentation system taken on a plane going through the center of the tip. 
The system is discretized to a lattice with L = 80. The tip and the indentation are taken to be half 
ellipsoids with vertical semi-axes a = 18 and horizontal semi-axes 6 = 15 and 6 = 19 respectively. 
In the figure / = 6. 

In the absence of a general method to compute dispersion forces in arbitrary geometries, 
two approximations are commonly used. The first is the proximity force approximation 
according to which the interaction between surfaces of any geometry is found by assuming 
that the surfaces are locally fiat. One then sums an effective, local free energy of the form 

over the surfaces where h is the local distance between bodies and A is the Hamaker constant. 

For the tip of Figure [1] we compared results obtained via proximity force approximation 
with our numerical results for ratios of the dielectric constants varying from = 1.05 to 

= 100. In order to construct the proximity force approximation, for each lattice site of 
one surface we find the closest lattice site on the other surface and build a minimal distance 



map. We then use the generalization of eq. ([6]) to periodic boundary conditions [l5| to 



evaluate the interaction at each distance in the map and add all contributions to obtain 



FIG. 2: Minimal distance map for / = 6 (top) and I = 20 (bottom). Black corresponds to 
the minimal distance, while white corresponds to the maximal distance. As the tip moves up 
the geometry of the closest approach changes from a central disk to a ring from the rim of the 
indentation. Mesh corresponds to the discretization lattice. 

the final free energy. The maps for the two chosen gaps are shown in Figure El We find 
that the free energies computed with the proximity force approximation can deviate to up 
to 40% from the full numerical evaluation. In particular the proximity force approximation 
performs the worse for high dielectric contrast (r^ = 50) and when the two surfaces are 
close together. Results for this system are shown in Figure [31 When the / = 20 proximity 
force approximation performs better with deviations from the full evaluation of 8% for large 
contrast. 

The second method, widely used in computer modeling, is an additivity assumption 
according to which the overall interaction is found by summing pairwise interactions of ele- 
mentary components. A typical example of such approximation performs an integration of 
the long-ranged tail of Lennard- Jones potentials between infinitesimal constituents of 



the interacting macroscopic objects 



22j. More generally, if one assumes a pairwise interac- 



tion between elements of a dielectric of the form V{r) = av{r) where a characterizes the 
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FIG. 3: Ratio between the free energies computed via proximity force approximation and our 
numerical results as a function of the dielectric contrast between the two materials for / = 6. 

strength of the potential, then the long-ranged part of the interaction energy between two 
macroscopic bodies can be written as 



where the coordinates R are the relative positions of the bodies and the function G encodes 
all the geometric information. 

We notice that both eq. ([6]) and eq. ([7]) display a factorization property (they are expressed 
as a product of material and of geometric terms) so that if we for instance consider the ratio 



for two different geometries we find a result independent of the material property. Since 
this factorization property is true in two very different limiting approximations it seems 
interesting to study the degree to which it remains valid over a wide range of dielectric 
contrasts in a geometry which is far from planar. 

Surprisingly, Figure HI we find that the ratio of free energies computed at different dis- 
tances exhibits only a moderate variation with the dielectric contrast. The robustness of this 
result has been tested by changing the system size and the accuracy of the tip/indentation 
discretization. Our study suggests a way of inferring free energy values at high dielectric 
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FIG. 4: Test of the factorization property, eq. ([8|). The ratio of interaction energies for two 
separations of the tip-indentation system of Figured] is plotted as a function of r^. The ratio varies 
by just 13% over the whole range of dielectric contrasts, despite the great variation in geometry 
displayed in Figure [21 



contrast once the results at low contrast are known. It is sufficient to determine the inter- 
action at low e to obtain the whole set of measurements for the desired system. Our results 
show that for classical, thermal interaction, the factorization method we propose, works 
better than proximity force approximation when large surface deformations are present and 
for high dielectric contrasts, both instances where proximity force approximation performs 
poorly. This approach can be combined with analytic approximations that compute free 
energies for arbitrary geometries in a small dielectric contrast expansion 14i\ . 



V. NON-RETARDED QUANTUM INTERACTIONS 

The quantum interaction between two materials is classified as non-retarded when the 
interaction is instantaneous, or retarded when one must take into account the finite propa- 
gation speed of electromagnetic radiation. In the first case the decay of interactions between 
two atoms is 1/r^; in the second the interaction falls as 1/r^. The characteristic crossover 
distance between these two forms is determined by the wavelength of typical spectral fea- 
tures that dominate dielectric response. This is often a feature in the ultraviolet, leading to 
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a cross-over length of tens of nanometers for most materials. 

We now generalize our treatment of the interaction to the short distance, non-retarded 
regime. The systems for which this effect is dominant are limited to the nano-scale, however 
we will show that this case is particularly simple to treat with the methods developed above. 
We leave the generalization to retarded interactions, which require the correct treatment of 
the propagating electromagnetic field, for a future publication. 



A. Quantization 

Since in the energy eq. ([1]) the field p represents a microscopic polarization vector we can 
study its dynamics by adding the kinetic energy. 

T = p(r)pV2 (9) 

where p{r) is a mass density. 

The thermodynamics of a quantum system are particularly simple to treat with the 



method of path integral quantization [23| . The potential and kinetic energies are combined 
in an effective statistical weight for an ensemble of N identical replicas of the original system; 
these replicas are coupled in the time/temperature direction by harmonic springs. The exact 
quantum partition function is then generated in the limit of large N. 
The effective action at each time slice, m is 

U^ = l dV + rU, (10) 

with r = (3/N. The harmonic coupling between the replicas comes from the presence of the 
terms in m and m + 1 in eq. (|TOl) arising from the kinetic energy. We now perform a Fourier 
transform in the time/temperature direction and find that the energy can be written in a 
form which is identical to eq. ([1]) if we define a frequency dependent susceptibility. 

^ = - + ifl(2-2cosa;) 

with Fourier frequency u = 27m/ N. 

We again perform the transformation from eq. ([1]) to eq. ([3]) by introducing an integral 
over the potential. After integrating over p the partition function of the quantum system is 
given by a product 

N-l 

Z=ll Z{e{n)) (11) 



n=0 
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where in the Z{e) of eq. (jlj), the dielectric function is replaced by 



= 1 + . , ....o o.Jo,^/.n^/,2,,2^2 (12) 



with uo = 1 / a/xoWpW- If introduce the Matsubara frequencies Un = 2m: / (3 we find 
that for large eq. (fT2i) simplifies to the single pole approximation 

often used ^\ to fit the dielectric properties of material in calculations of the Hamaker 
constant. Passage to the limit eq. f|T3l) requires that 

A^' > (^o/5)'e(0) (14) 

in order that the high frequency limit t{uj) — > 1 is correctly reproduced. 

This path integral formulation gives the full, combined thermal and quantum contribution 
to the interaction potential. It does, however, become less efficient at low temperatures; at 
exactly zero temperature the free energy should be evaluated by direct numerical quadrature 
over frequency. 

B. Implementation 

Each contribution to eq. ffTTj) at frequency n requires the evaluation of a single matrix 
determinant identical in form to that evaluated in the classical interaction. We measured the 
absolute discretizion error by studying the free energy of interaction between two parallel 
slabs for different values of N . We followed the procedure of [l^ using e(0) = 5, hw^jS = 40. 
If the largest contributions to the free energy comes from the ultra-violet range of the 
spectrum our choice of the parameters corresponds to a system near room temperature. 
Quantum effects then dominate the interaction, but thermal effects are correctly included 
in the evaluation. 

From eq. (fT4l) we require that ^ 90. We performed evaluations of the free energy for 
80 < < 260, Figure El Beyond A^ = 180 the free energy can be fitted by the Richardson 
form JF = J^oo + b/N'^ 25|, where J-'oo and b are fitting parameters. Using values of A^ greater 
than 180 we estimate JF^. We then evaluate the error generated for a specific value of A^. 
We find an agreement that varies from 93% for A^ = 80 to 99.0% for A^ = 260. We adopted 

12 



= 240 for the evaluation of the free energy landscape described below. For this value 
-^oo/^24o = 0.987. Higher accuracy and faster convergence (in are possible if we work 

with two values of N and extrapolate to J-'oo- 
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FIG. 5: Convergence of free energy differences with copy number A^. We evaluated the free energy 
difference of two parallel plates separated by ? = 2 and / = 10 for N varying from 80 to 260. 
The free energy is plotted as a function of N~'^. From values of above 180 we extrapolate 
Foo = 2.051 X 10-^. L=48. 

We then evaluated the free energy of a system composed of a sharp point over a surface 
with regular wells as a function of the horizontal position. Given that typical small force 
microscopy tip sizes are of the order of lOnm — 50nm forces measured experimentally by 
this technique will generally fall into the crossover between the retarded and non-retarded 
regimes. However, there do exist extreme cases of tips of radii of 2nm (super sharp silicon 
tips), where this case could be relevant, even if we are beginning to be close to the atomic 
scale where a continuum dielectric description is more difficult to justify. 

We scanned a region of 12 x 12 lattice points, covering the area of one well and its 
surroundings. Results of the free energy landscape are shown in Figure [61 The free energy of 
each point J-'{x, y) is evaluated for = 240, for which, due to the symmetry e{n) = e{N — n), 
only N/2 + 1 = 121 determinants are needed. Firstly, the matrices (12^ x 121 ~ 17,500) 
were built using Matlab, and stored on disk. This took approximately two days on a single 
workstation. Then, the free energy evaluation was performed on a cluster of nine processors. 
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and took approximately 2 days to complete, with the evaluation of the interaction at each 
position taking about 3 hours. The energy landscape we obtained reflects the underlying 
well profile, but with smoothed features. 




FIG. 6: Quantum free energy landscape for a system composed of a sharp, conical, tip near a 
surface with regular wells, side and depth / = 6, separated by 6 lattice spaces in both x and y 
direction. Box size L = 48. A portion of the system of dimensions 12 x 12 is scanned and the 
result reproduced periodically. Top surface: the free energy landscape at one lattice space above 
the top of the wells. Lower surface indicates the positions of the wells. 



VI. DISCUSSION 

Quite fine and useful discretization of fluctuating dielectrics can be studied using methods 
based on Cholesky factorization of a sparse determinant. One can study the full multibody 
dispersion interaction at nanometric length scales, a scale suitable for the study of macro- 
molecules and nano-particles, in non-trivial three dimensional geometries. Conventional 
methods for interpreting such systems use extensive modeling with force flelds that contain 
many simplifying assumptions which can now be checked (in the near fleld regime) against 
explicit numerical results. 



We now compare our approach with recent numerical work 
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ll| . The authors work with 



the full discretized Maxwell equations which enables them to consider the full retarded 
interaction between bodies; they calculate the stress tensor rather than the free energy that 
we chose to evaluate. The method requires the calculation of Green functions which are 
then integrated over a closed surface surrounding the body of interest. The authors propose 
several different methods of solving for the Green functions, including fast multipole and 
multigrid methods, which they argue can solve the problem in with memory S ^ V and time 
Te ~ ^/^-i/'^. They work with modest system sizes of 20 x 40 in simplified 2 + 1 dimensional 
geometries, suitable for studing grooved surfaces. 

In practice they used a conjugate method which requires a number of iterations which 
increase with the systems size as V^^'^ so that their actual implementation scales like 
our own as Tg ~ V"^. It could have been useful to compare real computing times, and 
not only theoretically derived asymptotic behaviors, but unfortunately the authors do not 
report these data. Any consideration of the choice of algorithm must take into account 
the prefactors in these laws; fast multipole methods have been abandoned in applications 
such as molecular dynamics due to the enormous prefactors in the (apparently favorable) 
asymptotic scaling. 

For system sizes comparable to those used in this paper one finds that the number of 
fioating point operations needed to perform the Cholesky factorization is comparable to 



^fiop = Q-5V'^ [27|. Multigrid methods while being asymptotically fast can require hundreds 
of iterations in order to converge 28|] . It thus seems possible that they perform less well than 
the methods of the present paper for moderate system sizes. It would be most interesting 
to study the crossover point in the efficiency of the various methods. 

When we restrict ourselves to systems in 2 + 1 dimensions the Cholesky factorization we 
use also improves in performance. Storage requirement in this case is only S ~ A^log(A^), 
and computation time is Te ~ N'^^^. The interaction of grooved surfaces is studied by 
Fourier transforming in the uniform direction before performing the factorization p^. In 
such a system evaluation of the interaction of a system discretized to a lattice of dimensions 
500^ can be performed in 30 minutes on a 3.2GHz 64-bit Xenon workstation. 

The main physical problems that are still not possible to study with this method have 
length scales in the range 20nm — 10/im where the retarded interaction dominates. This 
requires a different discretization strategy based on the full Maxwell equations and will be 
considered in an future paper [29 1. 
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